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ABSTRACT 

Aims. Some basic structures in planetary nebulae are modeled as self-organized magnetohydrodynamic (MHD) 
plasma configurations with radial flow. 

Methods. These configurations are described by time self-similar dynamics, where space and time dependences 
of each physical variable are in separable form. Axisymmetric toroidal MHD plasma configuration is solved 
under the gravitational field of a central star of mass M. 

Results. With an azimuthal magnetic field, this self-similar MHD model provides an equatorial structure in 
the form of an axisymmetric torus with nested and closed toroidal magnetic field lines. In the absence of 
an azimuthal magnetic field, this formulation models the basic features of bipolar planetary nebulae. The 
evolution function, which accounts for the time evolution of the system, has a bounded and an unbounded 
evolution track governed respectively by a negative and positive energy density constant H. 
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1. Introduction 
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Planetary nebulae are believed to be intermediate-mass stars that are expelling their outer layers, at slow 
velocity, during the asymptotic giant branch (AGB) phase. As the temperature of the star increases, the ejected 
material becomes tenuous but moves with increasing speed. The hydrodynamic interaction of these slow and fast 
winds, at the intersection of their trajectories, generates shock waves that create the characteristic morphologies 
of planetary nelulac (Dyson and de Vries 1972, Kwok et al 1978). The outward-shock compresses and heats 
the slow dense wind that is ahead, and the inward-shock decelerates and heats the fast wind behind. This 
interacting-wind model reproduces the spherical features of planetary nebulae. In contrast, the elliptical and 
bipolar features may be created by a dense toroidal cloud in the equatorial plane (Kahn and West 1985, 
Mcllcma et al 1991). Although toroidal clouds are imaged by high-resolution instruments, how such structures 
formed during the AGB phase is still a debated issue. Besides spherical, elliptical, and bipolar features, high- 
resolution images have revealed point-symmetric features (Miranda and Solf 1992, Lopez et al 1993, Balick et 
al 1993). These point-symmetric features, plus the detection of significant magnetic fields inside central stars 
(Jordan et al 2006), appear to call for a magnetohydrodynamic (MHD) approach to the study of planetary 
nebulae (Pascoli 1993, Chevalier and Luo 1994, Garcia-Segura 1997, Bogovalov and Tsinganos 1999, Matt et 
al 2000, Gardiner and Frank 2001). In this paper, we present such an MHD analysis, which is able to recreate 
axisymmetric structures such as an equatorial plasma torus and a bipolar nebula. We consider an MHD plasma 
that simulates a radial-flow explosion in spherical coordinates, and we solve for solutions that are self-similar in 
time (Sommerfeld 1950, Landau and Lifshitz 1978). 

Self-similar MHD analysis was pioneered by Low (1982a, b, 1984a, b), in astrophysical stellar-envelope appli- 
cations. The association of self-similar solutions to eruptive phenomena can be best demonstrated by the gener- 
ation of self-similar shocks in powerful atmospheric detonations. The shocks are the results of self-organization, 
through dissipative processes, following chaos. The shock is built up as the gas expands outward, and it takes its 
fully developed form beyond a certain radius. These ordered structures could have a scale invariance that lowers 
the dimensionality of the time-dependent fluid system, and could exist between one spatial coordinate and time. 
This invariance is identified as the self-similar parameter, and similarity is temporal. Time evolution would be 
a separable factor, and, as a result, the time and space dependencies of self-similar variables are represented 
in separable form. For steady-state systems, the scale invariance could be between one spatial coordinate and 
another, and the similarity is said to be spatial. 

Self-similarity is often regarded as a method to derive specific types of time-dependent or steady-state so- 
lutions. In this paper, we consider self-similarity to be a property of self-organized states, at least some of 
which are derived from turbulence. It is essential to note that the appearance of self-organized states, in fluids 
and magnctofluids, owes itself to the existence of multiple quadratic invariants in the absence of dissipations 
(Hasegawa 1985, Biskamp 1993). In the case of MHD fluids, the invariants are the total energy density (mag- 
netic, plasma, and thermal), magnetic hclicity, and in the case of incompressible fluids, cross-helicity. In the 
presence of dissipation, these invariants could undergo constrained changes, for example, from an MHD con- 



Tsui: Equatorial Plasma Torus 3 

figuration, permitted by the MHD equations, to another topologically accessible configuration of lower energy. 
It is important to point out that these configurations are isolated in ideal cases. The continuous topological 
transformation from one to another configuration is only possible through dissipative paths such as magnetic 
reconnections. Mathematically speaking, this evolution amounts to the application of the variational principle, 
on the total energy density, under the constraint of global magnetic hclicity conservation (Hascgawa 1985). 
The exact processes of dissipation, in taking the system to self-organization, need not be specified under the 
variational description. For this reason, MHD systems have the tendency to develop self-organized states with 
structural stability. 

Since the self-similar method solves for the end configurations, directly from the governing equations, there 
is no need to specify the initial conditions. The plasma would find its way through magnetic reconnection, 
and other dissipative means, to reach one of the possible end configurations. The only condition is that the 
initial configuration is dissipatively transformable to the final configuration. If the initial conditions are close to 
those of the end configuration, the time to reach the end configuration is brief. Naturally, there are many end 
configurations, and only one of them will be selected by dissipative reorganization. For this reason, we cannot 
predict which self-similar state the system will reach. We can only match a system to a given self-similar solution 
that bears some resemblance. 

One of the fundamental obstacles to visualizing self-similar solutions is an attempt to associate a given 
initial configuration with a specific self-similar end solution. It is possible to predict the final state of ideal 
MHD systems, in particular, by iterating equations of motion for an initial configuration. In the presence of 
dissipations, it can, however, be difficult to predict the final state, because of intermediate, chaotic events that 
occur while states are being reorganized. The initial configuration could undergo topological changes, through 
dissipative processes, such as magnetic reconnections, subject to the constraint of quadratic invariants. To 
describe the self-organized states arising from chaos, we have to override the forward time iteration process and 
attempt to identify the end configurations directly from the governing equations of the system. Although these 
end configurations have originated from dissipative reorganization through chaos, it is important to emphasize 
that these end configurations are determined by the ideal governing equations. Ideal MHD cannot, however, 
determine the dissipative path taken to reach one configuration from another. As a consequence, we are unable to 
recover the initial conditions by following the self-similar configurations backward in time. The initial conditions 
are lost from memory. 

Recently, time self-similar analysis was applied to represent extragalactic jets as a polar-launched ejection 
event (Tsui and Serbeto 2007), which differs from the classical accretion-ejection space self-similar steady state 
MHD transport model (Blandford and Payne 1982). A self-similar description has been applied to different 
problems in physics, such as atmospheric ball lightnings (Tsui 2006, Tsui ct al 2006), and interplanetary magnetic 
ropes (Osherovich et al 1993, 1995, Tsui and Tavares 2005). Here, we follow the methods of time self-similar 
MHD analysis to model plasma torus, by developing solutions that are launched onto the equatorial plane. 
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2. Axisymmetric Self-Similar MHD Scalings 

The basic MHD equations, in Eulerian fluid description, are given by 

|? + v-g™) = o, (l) 

r dv , „. , _ GM 

P{ ~dt + ^ ' )V} = J x B Vp — p—r , (2) 

B Ti 

¥ = -Vx£ = Vx(„xB), (3) 

VxB = f iJ , (4) 

V • .B = , (5) 

■!■(—) + («• V)(— ) = , (6) 

where p is the mass density, v is the bulk velocity, J is the current density, B is the magnetic field, p is the 
plasma pressure, \i is the free space permeability, 7 is the polytropic index, and M is the central mass, which 
provides the gravitational field. 

We consider a radially expanding plasma, and look for self-similar solutions in time, where time evolution 
is described by the dimensionless evolution function y(t). For this purpose, it is most convenient to think of 
Lagrangian fluid description, and consider the position vector of a given laminar flow fluid element r(t). Under 
self-similarity, the radial profile is time invariant in terms of the radial label 77 — r(t)/y(t), which has the 
dimension of r. Physically, 77 is the Lagrangian radial position of a fixed fluid element. With a finite plasma, 
the domain of 77 is bounded by mass conservation 

< Vint < V < Vext ■ 



As for the plasma velocity, we consider self-similar structures deriving from a spherically symmetric radial 
velocity that can be written as 

{-1.0,0} -{..0.0}. (7) 



Our self-similar parameter 77, defined through the Lagrangian fluid label, explicitly represents the fluid velocity 
by the time evolution function y(t). This evolution function will be solved self-consistently with respect to the 
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spatial structures of the plasma. We emphasize that self-similarity, as a method, can be applied in different 
ways other than the one we use here. For example, Lou and his collaborators have treated an aggregating 
fluid under its self-gravitational field with a similarity variable x = r(t)/at, where a is the sound speed, for 
isothermal fluid, 7 = 1, (Lou and Shcn 2004, Bian and Lou 2005) and for a polytropic gas, 7 > 1, (Lou and 
Wang 2006, Lou and Gao 2006) to study relevant astrophysical phenomena. Extensions to a magnetofluid have 
been considered by Yu and Lou (2005) and by Lou and Wang (2007). Because of the linear dependence on 
time, this similarity variable x refers to a reference frame moving at speed a, which is different from the radial 
plasma flow velocity v. Furthermore, different from our similarity variable ij, x here is not the Lagrangian label 
of a given fluid element. For this reason, the convective derivative remains explicit in the x representation. As 
a result, the similarity variable of Lou amounts to finding the plasma structures in an adequate moving frame 
in the Eulerian x fluid description. This resembles the analytic technique of going to a moving frame to look 
for stationary profile solutions for nonlinear phenomena such as nonlinear Alfven waves, solitons, etc. Because 
of this fundamentally different definition of self-similarity, the nature of the phenomena intended to describe 
is different. Our Lagrangian label self-similarity parameter is aimed to find spatial plasma configurations, and 
to determine the radial plasma flow velocity, consistent to the spatial configurations, through the evolution 
function. 

Since we are considering an isotropic radial plasma flow, a natural solution would be a hydrodynamic 
one-dimensional expanding plasma, with radially dependent mass density p and plasma pressure p, and with 
B = and J = 0. Nevertheless, this one-dimensional hydrodynamic solution is highly unlikely because 
magnetic field fluctuations can be generated from current density fluctuations, even in the absence of a 
pre-existing magnetic field. With the magnetic fields, which are basically a two- or three-dimensional structure, 
coupling to the plasma will generate likewise two- or three-dimensional p and p. To examine the possible spatial 
structures of MHD plasmas, we consider a two-dimensional case with azimuthal symmetry in <f>. In this case, 
the magnetic field, through the vector potential A, can be expressed as 

This enables the magnetic field, the current density, and the electric field, to be expressed as two scalar 
functions P and Q respectively 

B = -^—A+-i- a (P), ~S-(P),+Q} = VP x V0 + QV</> , (8a) 
r sin 6 r 06 Or 

1 r 18Q _8Q _8^P_l_ . a 8 , 1 dP^ 
^ J ~ rsin9 i+ r 89' 8r ' 5r 2 r 2 Sm 89 [ sm9 89 j) ' ( j 

1 8P 

E = -(v x B) = —^{0,vQ,v— } . (8c) 
r sm 9 or 
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The independent variables are transformed from (r, 6, t) space to (77, 0, y) space. We determine the explicit 
dependence of y on each one of the physical variables with this radial velocity using functional analysis. First, 
making use of Eq.(7), Eq.(l) renders 

dp .dp dp, ,dv 2v, dpdy 3p dy .dp 3p.dy „ .„ . 

i + ^8-r^ = ti+v£ )+p{ - + -) = ^| + f| = (^ + f)J=0. (9a) 

To reach the second equality, we note that the first bracket, in the first equality, corresponds to the total time 
derivative of an Eulerian fluid element, which amounts to the time derivative of a Lagrangian fluid element. As 
for the second bracket, it can be reduced by using v = dr/dt = rjdy/dt and dv/dr = (l/y)(dy/dt). Solving this 
equation for y scaling, by separating the time part, gives 

p(r,t) = ^p( V ,9) . (96) 



As for Eq.(6), with aoF = (p/p 1 ) where «o is a constant that carries the physical dimension so that F is a 
dimensionless function, it follows 

dF dF dFdy , . 



(-£) = a F(r,t) = ^aoF^e) . (106) 



As for Eq.(3), using the representation of Eq.(8a), the magnetic functions P and Q are 

dP dP dPdy . 

P(r,t) = ^P(v,e) , (116) 



dt dr W > dy dt y dt [ dy y ' dt 



(12a) 



Q{r,t) = ^Q(V,0) 



(126) 
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3. Self-Similar Formulation 

We have reduced the general set of time-dependent ideal MHD equations, Eqs.(l-6), to a set of self-similar 
equations with appropriate time scalings, Eqs.(7-12), and have sought solutions of Eq.(2) and the time evolution 
function y(t). The general ideal MHD set of equations has nonlinear terms of the convective type (v ■ V) in 
Eq.(l) and Eq.(2), and interaction type (J x B) in Eq.(2) and (v x B) in Eq.(3). By using the fluid label 
description, the (v ■ V) convective terms are absorbed into the Lagrangian time derivative representation. The 
(v x B) interaction term is also absorbed, through an adequate representation of the magnetic field, by scalar 
functions P and Q. The structure of the nonlinear terms, absorbed into the Lagrangian fluid label formulation, 
will appear in the r\ profile of the system. The remaining task is to solve Eq.(2) for the self-similar r\ structure, 
under the presence of the (J x B) interaction term. 

We note that the equation of a$F = {p/ p 1 ) and the equation of P are of the same form, and conclude that 
F = F(P) is a functional of P only, or F — F(P). The fact that two equations have the same form does not 
automatically imply that the solution of one equation is a functional of the solution to the other. We can make 
this assertion only when considering self-similar solutions. We proceed to write the (r], 9) dependences in terms 
of (?7, P(r], 6)) in p and p, to obtain 



p = p( V ,P) , (13a) 
p = ±. ao F(ri,0)p'{ V ,e) = ^a Q F(P)p^f 1 , P) - ^p(v,P) • (136) 

Furthermore, since (p/ ' p 1 ) — aoF(P) is a function of P only, the rj and P dependences should be in a separable 
form in both 



P(V,P) = PoPi(v)P2(P) , 



p(r],P) = PoPi(v)P2(P) ■ 



The dimensions of mass density and pressure appear explicitly in p and po respectively. In this form, pi(r]), 
pi (P) , pi (rj) , and p2 (P) are dimensionless functions. By taking oto = po/ Pq with the physical dimension, we have 

Pi fa) = PKv) , (14a) 



p 2 (P) = F{P)pI{P) , (146) 
such that F should be a functional of P only. 
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Making use of Eq.(4) to eliminate the current density in Eq.(2), we derive the momentum equation, which 
has three components. The <j> component, which contains only the magnetic force, gives 



dPdQ_dP_dQ J_r^PdQ _ dP_dQ, _ 

~dr~~dJ ~ ~d6~d7 ~ ^ { ~£hj~de ~ ~de 7v ~ ' ( ' 

Q = Q(P) ■ (16) 

The function P plus the functional form of Q, therefore, determines the magnetic fields. As for the 9 component, 
with p = p(r), P)/y Sl , it reads 

S + J 2 sin 4^S) + «i + ^ 4 " 3 v-^n = o . (it) 

orf r\ z oO smfl oO dP oP 

We remark that the first three terms of this equation represent the nonlinear force-free field equation with 
Q — Q(P). In the particular case of a linear functional Q — aP, this equation describes the linear force-free 
fields, which can be verified easily from Eqs.(8) with fiJ = aB. Such a force-free magnetic configuration is the 
marker of self-organized plasmas, and is the result of the variational principle that minimizes the energy density 
of the magnetic field, under the constraint of global magnetic helicity conservation. This justifies the use of 
self-similarity method to describe self-organized configurations. The last term of Eq.(17) is the plasma pressure 
term. For an isotropic pressure independent of P, this last term would be null and such plasma pressure p(r]) 
has no effect on the magnetic field which remains force- free. For an axisymmetric pressure p(r],P), the above 
equation would be independent of the time evolution function y should 7 = 4/3, and it would describe a plasma 
pressure balanced axisymmetric magnetic field. Such a pressure-balanced self-organized configuration could be 
derived using the variational principle, by minimizing the energy density that contains not only magnetic energy 
but also plasma thermal energy under the same constraint. The solution of this equation of P, with different 
representations of Q = Q(P) and p = p(P), is the core problem of self-similar description. 
As for the r component of the momentum equation, it reads 

r dv dv, d . n . ... GM 



dt dr dr ' r 2 

I, 1 , 2 dP t d 2 P 1 . n d , 1 dP. ^dQ^ dP d . „. 
yu rsmW or or z r z oO smu 08 oP or OP 

The term dp/dr on the left side refers to the total radial derivative of the plasma pressure, which includes both 
an explicit and implicit dependence in P. The right side of the radial equation provides an expression for the 
magnetic force. Making use of the meridian component, Eq.(17), the right side is equal to the implicit part of 
the radial pressure gradient, which cancels the same term on the left side. This leaves only the explicit radial 



pressure gradient 
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,dv dv, d , „, GM d 2 r d GM , <nI , 



In terms of self-similar parameters, the radial equation reads 

,.£ + V*4<* + «£ } _0. (19) 
ai^ par] rj z 

With a as the separation constant, a 2 as a dimcnsionless constant, and 7 = 4/3, we derive 

Po 1/3/ n /I 9 GM X a , , 2GM a , , „ q. ,„„ , 

4^%) = (- V + _)=_(,» + -_)= -(,3 + 2,3) f (20a) 

p 2 (P) = a 2 p 2 (P) , (206) 



The first of these three equations provides the mass density profile of the radial label. The second one specifics 
the P dependence between mass density and plasma pressure. The third one defines the time history of the 
evolution function. 
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4. Axisymmetric Equatorial Structures 

After describing our self-similar formulation, we solve Eqs.(14), (17), and (20), for axisymmetric plasma 
configurations. We begin with comparing Eqs.(14b) and (20b), and conclude that 

F(P) = a 2 p\-<{P) . (21) 

The solution to Eq.(17), for the magnetic field function P, depends on the functionals Q and p2- This is where 
different types of self-similar configurations can be produced. Although Low pioneered the use of self-similar 
MHD in astrophysics, he did not explore the full potential of the method. We remark that self-similar solutions 
are highly dependent on the functional representations of Q and p2- Under a given geometry, there are different 
self-similar solutions, with different characteristics to represent different self-organized plasmas, depending 
on how the magnetic fields and plasma parameters are modeled. To study configurations on and about the 
equatorial plane, we consider 

Q = aP , (22a) 



MP) = (P + C) > • ( 22b ) 

For example, if we assume a quadratic function to have 

p 2 (P) = (P 2 + C) > , 

we derive a polar ejection configuration that tends to collimate magnetic fields and plasma density, onto the 
polar axis, as plasma pressure builds up (Tsui and Serbeto 2007). Using these functional representations in 
Eq.(17), we derive 

2 d 2 P . J, 1 SP. , , 2 - 
or) 2 oO smfc> of) 

= -/ipo?7 4 Pi(?7)sin 2 6 = -UPo^Pi^iv) si n 2 # • (23) 

To solve this equation, we separate the variables by writing P(j],9) = R(rf)Q(6). Assuming x = cos9 and 
n(n + 1) to be the separation constant, we have 

(1 _ X 2)^W + n{n + 1)e{x) = } (24fl) 
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^0 + [M 2 - »(» + W = -^V 4 pt%)^ ■ (246) 

The first equation provides the well-known solution, in terms of the Legendre polynomial P n (x) , 

Q(x) = C a (l-x 2 )^-^ = -n(n + l)C f ' P n {x)dx . (25) 
ax h 

We note that the second equation is not in a separable form. Nevertheless, should we choose n = 1 with 
Pi(x) — x, such that Q(x) = (1 — x 2 ), where we assume that Co = 1, Eq.(24b) is separable. The general 
solution is given by a homogeneous solution Rq (rj) and a particular solution R\ (77) 

R(ri) = Ro(v) + Mv) = 4oh)jnH + iJiW , (26) 

where jn(arj) with n = 1 is the spherical Bessel function regular at r\ — 0. 

The coefficient A of the homogeneous solution carries the amplitude and dimension to reflect the magnetic 
field through Eq.(8a). Since plasma pressure and mass density are expressed in terms of P viap2(-P) and f>2(P), 
we could normalize the entire self-similar solution, with respect to the magnetic field, by taking Aq = 1. With 
Pi(r]) given by Eq.(20a), denoting z — ar], and k 3 = 2(ary*) 3 = 2z 3 , the particular solution is described by 

z 2 ^ + [z 2 -n{n+l)]R 1 = -^(.^ J-£°)V + « 3 ) 4 = -A(z 3 + « 3 ) 4 . 
dz z a A %a 2 po 

We note that the coefficient A, which has two factors, is dimensionless. The first factor (fipo/a 4 ) ~ /3/z 2 is 
related to the plasma /3, which is the ratio of plasma pressure to magnetic pressure. The second dimensionless 
factor with GM is basically the ratio of the gravitational potential energy density to the plasma pressure. 
The right side of this equation can be expanded into binomial terms, and can be solved using power series, 
by taking one term at a time. The particular solution is, therefore, a superposition of five sub-solutions 
#1(12)1 #1(9)1 #1(6)) #1(3) > an d #1(0)) where the bracketed number in the subscript denotes the power of z 
in the binomial term, on the right side. We denote the coefficients A^ 12 ^ = A, A^ = 4k 3 A, A^ = 6k 6 A, 
^4(3) — 4 K 9 A, A(°) = n 12 A. We solve each binomial term of the particular solution using a power series of the form 

Using standard power-series techniques, we derive solutions of the form 



#1(12) ( z ) = a 2Z 2 + a^z 4 + a e z 6 + a 8 z 8 + a w z 



(27a) 



12 
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where a w = -A (12 \ a$ = -(10 x 9 - 2)a w , a 6 = -(8 x 7 - 2)a 8 , a 4 = -(6 x 5 - 2)a 6 , a 2 = -(4 x 3 - 2)a 4 , 
and ao = -(2 x 1 - 2)a 2 = 0, 

i? 1(9) (z) = a 9 z 9 + anz 11 + a 13 z 13 + a 15 z 15 + , (276) 

where a 9 = -A< 9 )/(9 x 8 - 2), a 1± = -a 9 /(ll x 10 - 2), 013 = -On/(13 x 12 - 2), a 15 = -oi 3 /(15 x 14 - 2), 
a m = -a m - 2 /(m x (to - 1) - 2), 

#i(6)( z ) = fl22 2 +a 4 z 4 , (27c) 

where a 4 = -A^ 6 ', a 2 = -(4 x 3 - 2)a 4 , and a = — (2 x 1 - 2)a 2 = 0, 

#1(3) (2;) = a 3 z 3 + a 5 z 5 + a 7 z 7 + a 9 z 9 + , (27d) 

where a 3 = -A^/('S x 3 - 2), a 5 = -a 3 /(5 x 4 - 2), a 7 = -o 5 /(7 x 6 - 2), a 9 = -a 7 /(9 x 8 - 2), 
a m = — a m _ 2 /(m x (m — 1) — 2). As for i? 1 ( )(z), this solution has to be solved with series of negative powers, 
which describes the divergent nature of the gravitational field GM/rj at r\ = 0, although this limit is irrelevant 
here. We derive 

where a = 0, a 2 = a 4 = -(2 x 3 - 2)a 2 /(4 x 5), a 6 = -(4 x 5 - 2)a 4 /(6 x 7), a 8 = -(6 x 7-2)a 6 /(8 x 9), 

a m = — ((m — 2) x (m — 1) — 2)a m _ 2 /(TO x (to + 1)). We note that the homogeneous solution corresponds to 
force-free magnetic field configurations. Due to the plasma pressure in Eq.(17), there are particular solutions 
that maintain the pressure balance. 
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5. Magnetic Fields and Evolution Function 

After solving for the spatial structure of P, the magnetic fields are given by 

1 - , „ N 1 1 1 dP 11, ,de(x) , . 

y z y z rj sin 9 r\ 69 y z r\ z dx 

Be = - 2 Be{v,0) = -~ 2 —^ = \u/2 ^) > ( 28& ) 

y z y z rjsmv ar\ y z r\ drj (l-x z ) L / z 

= \B4 V ,9) - +L-L-Q = +^R(rj) — j— -y^0(x) . (28c) 

yz y z -qsmO y z f] {l—x z ) l l z 

The self-similar homogeneous solution i?o(?y), provided by Eq.(26), enables oscillations in r\ to be modeled, 
and vanishes at the zeros of the spherical Bcssel function. To understand the magnetic structures, we first 
set aside the particular solution Ri(r/). At the first zero of Ro(r]), we have B^{rf) = 0, B r (rj) = 0, and the 
only nonvanishing field is B${rf) = according to Eqs.(28). The meridian self-similar solution &(x), given by 
Eq.(25), oscillates in x. Together they describe the magnetic fields of Eqs.(28). Within this region of (rj, x), the 
topological center, defined by dR (r])/di] — and dQ(x)/dx — 0, has B r = 0, Bg = 0, and the only nonvan- 
ishing field is B$(rj) — 0. This is the magnetic axis. The field lines about the magnetic axis are, therefore, given by 

B r Bg Ba. 



dr rd6 r sin 9d(f> 

By axisymmetry, the third group is decoupled from the first two groups. For the field lines on an (r — 9) plane, 
we consider the first equality between B r and Bg, which gives 

P{r],x) = Ro(v)(l-x 2 ) = C . (29a) 



The nested field lines are given by the contours of P(r], x) on the (r — 9) plane. At the topological center, the 
magnetic axis, we have <d(x) at its maximum value, and i?o( ? 7) at its maximum, such that P(r],x) takes its 
maximum value. Since 

27rrsin6»B0 = 27raP , (2%) 



by Eq.(28c), the line integral of B^ on the magnetic axis is a maximum. We, therefore, have a sequence of 
plasma tori, due to the periodic nature of the spherical Bessel function. The presence of the particular solution 
Ri(rj) modifies this structure. Due to the divergent nature of Ri(rj), as we will see in the next section, only the 
first torus prevails. 
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As for the temporal part, the evolution function in Eq.(20c) is described by 

(f ) = [2(H + . (306) 

Here, H is an integration constant that is independent of time. To understand the meaning of a, we note that, 
using Eq.(20c), plasma acceleration in Lagrangian coordinates is 

dv d 2 y arj 

Tt =7] W = -^ ■ (31fl) 



A positive a corresponds to an outward, decelerating flow. The deceleration becomes smaller as y, or as r, 
becomes larger, and a equals the intensity of the deceleration, for a given radial label r\. To derive an expression 
for H, we multiply Eq.(30a) by prj 2 , and use Eq.(31a) to obtain 

2 1 2 , dv N 1 o , d , . GM N , , N 

pry 2 !/ = -pv 2 + (p— )r = -p« 2 - (^p(r, P) + p— )r . (316) 



The first term on the right side is the kinetic energy density of the fluid element, and the second term is its 
work on expansion, due to the explicit part of the pressure gradient. We note that the magnetic energy does not 
appear in this expression because it is cancelled by the work due to the implicit pressure gradient, as discussed 
in Eq.(18a). It is clear that H measures the total energy of the fluid element. 
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6. Force-Free Equatorial Plasma Torus 

We consider the explicit radial dependent part of the mass density, in Eq.(20a). Since p\(rj) is positive, a 
should also be positive. Differentiating with respect to 77 indicates that pi(rj) has a minimum at 

, GM , N 

ri = — . (32) 



Because we are considering a finite plasma, the range of the radial label 77 is finite. We assume that 77* is enclosed 
in the interval < rji nt < 77* < rj ext . For 77 > 77*, p~i(r)) increases because of the arj 2 /2 term, where a is defined 
by Eq.(20c) and measures the inertia of the plasma. A positive a corresponds to a decelerating plasma. Beyond 
the external bound is interstellar space, where plasma mass density decreases abruptly. In the presence of the 
gravitational term, GM/77, in Eq.(20a), the mass density p~i(r}) rises again for 77 < 77*. This gravitational term 
is singular at 77 = 0. However this singularity is not included, because the internal bound of 77 is at the surface 
of the star. 

To consider the function P(rj,x) — R(i])$(x), the homogeneous solution Ro{rj) = A (ar])ji(ar]), with 
A = 1, is shown in Fig.l, with the first node at z = (077) = 4.5. In the absence of plasma pressure, Ro{v) 
would be the only term for i?(?7). With z ext = arj ext = 4.5, this would correspond to the equatorial plasma 
torus force-free magnetic fields. As for the particular solution Ri(rj), we note that i?i(i2) (2), i?i(g)(z), Ri(q)(z), 
and i?!( 3 )(z) diverge as z goes to infinity, whereas Ri^^z) diverges as z goes to zero. These two extremes are 
naturally excluded, because z — is the center of the central star, and z — 00 requires an infinite plasma, which 
is unphysical. An adequate domain could be provided by the homogeneous solution plotted in Fig.l, which 
shows an axisymmetric equatorial plasma torus with < z < 4.5, if the particular solution does not significantly 
alter this structure. The particular binomial solutions depend on the parameter k 3 = 2GMa 3 /a = 2(ar]*) 3 , 
where a is a measure of the azimuthal magnetic field given by Eq.(22a). We assume the minimum z* — a?7» = 2, 
such that k 3 — 16. We note that the two factors (brackets) in 

A = /MPOw 1 1 PqGM 1 4 
( a i,[ 8a 2 z2 77, ?V 

in the equation of Ri(r]) are dimensionless quantities. The first factor (/ipo/ fl4 ) — P/ z * — 10 _1 can be expressed 
in terms of plasma /?, which is the ratio of the plasma pressure to the magnetic pressure. This factor is estimated 
to be 10 _1 . The second factor is basically the ratio of the gravitational potential energy density poGM/r)* to the 
plasma pressure po, plus other dimensionless multipliers to the fourth power. We believe that the gravitational 
potential energy density and the plasma pressure are approximately similar. We could take a 2 = 1, and the 
remaining terms would provide the factor 10~ 6 . Multiplying the two factors together gives A = 10~ 7 . 

Of the five parts of the particular solution Ri, R\n2) is unaffected by the choice of k 3 , and is plotted in Fig. 2. 
Other binomial parts Rug), Ri(q), Ri(3), plotted in Figs. 3-5 respectively, carry mixed contributions of plasma 
pressure and gravitational field, and are negligible within the domain. The last part i?i(o) in Fig. 6 is purely 
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gravitational. The neighborhood of the center of the star, z = 0, is excluded. These binomial parts depend on 
the choice of k 3 , but would have a negligible effect on R , within the domain of interest. The homogeneous 
and the particular solutions provide the complete solution for R(t]). Together with Q(x), they determine the 
function P(r], x), which governs the entire axisymmetric self-similar profile in magnetic fields, plasma pressure, 
and mass density. 

To consider the magnetic field profiles, we recall that the field lines on an (r — 0) plane are given by Eq.(29a), 
with R{z) = Rq{z) + R\{z). The presence of a plasma gives a diamagnetic effect, as shown by R\(i2){ z ) m Fig. 2, 
carrying an opposite sign with respect to Rq{z) of Fig.l. Examining Fig. 2, we note that i?i(i2) vanishes at 
z = 4.5 and z = 5.2. The first node of -Ri(i2) happens to be at the same location as the first node of R , also at 
z = 4.5. Furthermore, Figs. (3-6) show that i?i(g), Ri(e), Ri(3), Ri(o) are all negligibly small. Consequently, the 
general solution R{z) = Rq(z) + R\(z) = at z = 4.5. By Eqs.(28), we conclude that the radial and azimuthal 
magnetic fields, B r (z) and B<p(z), vanish at this location, while the meridian magnetic field, Bg(z), is maximum. 
Consequently, the domain < z < 4.5 constitutes a region for a magnetic torus with closed field lines. Since 
plasma is tied to the magnetic field lines, it circulates along and is confined within this torus. We note that the 
sign of R\ is negative, which indicates the diamagnetic nature of plasma. We note that i?i(i2) has a maximum 
at the same location as Rq. As a matter of fact, the functional form of i?i(i2) in the interval (0,4.5) is much 
the same as that of the force- free R$. Consequently, the general solution R is practically force- free, even when 
plasma pressure is taken into account, in Eq.(23). 

For simplicity, we neglect the diamagnetic effect, and take R{z) = Ro(z) to highlight the poloidal field lines 
in Fig. 7, which is described by first equality of the field line equation, with solutions in Eq.(29a). The addition 
of i? 1 ( 12 )(z) would change the contours slightly, but not the topology. Apart from these poloidal fields, there 
are associated toroidal fields given by Eq.(28c). The maxima of R and 6, dP/drj — and dP/d9 = 0, give the 
maximum of this toroidal field where B r = and Be = respectively. This corresponds to the magnetic axis of 
the plasma torus. Together, the field lines thread out nested magnetic surfaces enclosing the magnetic axis, one 
within another. The plasma pressure and mass density of this plasma torus have R(z) > 0. With Q(x) > 0, this 
warrants P(z, x) > to assure positive pressure and mass density. With C — in Eq.(22b), the plasma pressure 
and mass density vanish at z = 4.5 since R(z) = 0. It is important to note that, for z > 4.5, the particular 
solution R\(z) dominates the homogeneous solution Rq(z) such that the general solution R(z) never vanishes. 
The radial component of the magnetic field is nonzero. For this reason, the field lines are open to outer space, 
and the plasma is no longer confined. This outer open plasma structure is not shown in Fig. 7. 

For the evolution function, it is most convenient to write Eq.(30b) in the normalized form 




dr y a 



(33) 



with r = (2a) 1 ' 2 t as the normalized time. We assume that a is positive, and take y(0) = 1 to be the initial value, 
such that the Lagrangian label rj corresponds to the initial position of the Eulcrian fluid clement. Equation (33) 
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indicates an asymptotically convergent solution, with H = —\H\ < negative, which is shown in the lower curve 
of Fig. 8 with \H\/a = 0.5. This track represents an asymptotically stable radius for the equatorial plasma torus. 
On the other hand, if the energy density H = +\H > is positive, there is an ever expanding evolution track 
with a terminal velocity dy/dr = {H/a) 1 / 2 . This is shown in the upper curve of Fig. 8, which is inappropriate 
for AGB slow winds. 

To better understand the amplitude and the sign of H, we multiply Eq.(30a) by 2pr\dr\, and integrate to derive 

2H J m dn= (|) 2 / W ^-2 . (34) 

The first integral on the right side is positive definite, and corresponds to the kinetic energy of the expanding 
structure. Inside the second integral, there is a plasma pressure term, and a gravitational term. The gravitational 
term is always positive. Within the plasma pressure term, we have dpi(ji)/dr] < for 77 < 77*, since pi(ri) has 
a similar profile as pi{rj), according to Eq.(14a). This part has a forward plasma pressure force. For 77 > 77,, 
we have dpi(r))/dri > 0, and the corresponding plasma pressure force is backward. It is clear that the explicit 
part of the pressure profile pi(r/), which is described by the explicit part of the mass density profile p\{vi) in 
Eq.(14a), contributes to the sign of H. Rewriting Eq.(20a) as 

Pi (V) = -: — ) +(— )] . (35) 

4a 2 Po 2 77* 77 

and assuming that the coefficient on the right side has a value of unity, these two profiles are shown in Fig. 9, 
where we have taken r\ — 377* or z = 3z» = 6 as the upper bound in Fig. 9, as an example. The pressure profile 
in the positive gradient region r\ > 77* helps to cancel the effects of the negative gradient, over the region r\ < 77,, 
in the above equation. Assuming that the slow wind velocity is below the escape velocity, the first integral on 
the right side is less than the gravitational term of the second integral. This implies that H is negative, due to 
the overall negative sign on the second integral. The plasma torus is, therefore, asymptotically stationary. 

When deriving the full radial profile of the plasma pressure, we caution that Fig. 9 is only the explicit part, 
pi{rj), of the radial profile relevant to Eq.(34). Another radial contribution comes from the implicit part, p2(P)- 
The full radial pressure profile is a combination of the explicit and the implicit part, p(r],P) = popi(i])p2(P), 
where p 2 (P) = P(r],x) = R(r])9(x). Since pi(jy) = pj(r]) by Eq.(14a) and p 2 {P) = a 2 p2(P) by Eq.(20b), we 
could plot the full mass density profile, which reflects the plasma pressure profile with C = 0, as 

p( V ,P) = PoPi(v)P2(P) = Po[^- P -^?[\(-) 2 + (^)?-R(v)e(x) . (36) 

4«2 Po a z A z* z a 2 

With z* = 2 and considering i?(?y) = Rq(t]), the full radial dependence p~i(rj)Ro(r]), with the value of the 
coefficient in Eq.(36) taken to be unity, is shown in Fig. 10 for the plasma torus with a very pronounced 
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maximum at z = 3.7, and zero at z = 4.5. The corresponding mass density spatial contours of p(r],P) are 
shown in Fig.ll. 



Tsui: Equatorial Plasma Torus 19 

7. Bipolar Planetary Nebulae 

In addition to the equatorial plasma torus, our self-similar MHD formulation is able to model planetary nebulae 
of the bipolar type. For this purpose, we consider 

Q = , (37) 

to warrant = 0. Following the same analysis, the radial function becomes 

R(r)) = Ro{r)) + Ri(r]) = A^ 1 - a 12 z 12 - a 9 z 9 - a 6 z 6 - a 3 z 3 + a a . (38) 

The coefficients are (12 x 11 - 2)a 12 = A (12) = A, (9 x 8 - 2)a 9 = A^ = 4k 3 A, (6x5- 2)a 6 = A^ = 6k 6 A, 
(3 x 2 — 2)a3 = A^ = 4k 9 A, and 2a = A^ = n 12 A. Although we have continued to use z = arj as the 
normalized radial label, the normalization factor a here has no connection to Q because we now assume that 
(5 = 0. We could imagine that a = k, such that z = krj. We notice that the particular solution consists entirely 
of negative terms, apart from the final term, whereas the homogeneous solution is positive. With Aq = 1, 
A = 10 , z* — 2, k 3 — 2z 3 = 16, the general solution R(i]) shown in Fig. 12 becomes negative at z = zq = 4.55. 
This value of Zq happens to be close to the roots of Ro(z) and R\{z), of the equatorial plasma torus case, 
where Q = aP ^ 0. This fact is purely a coincidence. The full radial profile of the plasma density, described 
by Eq.(36), explicit and implicit parts, is shown in Fig. 13. Since plasma density is positive definite, z < zq 
defines the physical domain of our self-similar solution. The bipolar field contours are shown in Fig. 14. The 
radial component vanishes at zo, B t (zq) = 0, which is the boundary of the expanding bipolar structure. 

To conclude, we have represented the AGB wind by an exploding ideal MHD plasma. The explosion is 
modelled by an isotropic radial plasma velocity in spherical coordinates. Through dissipation, this exploding 
plasma is believed to undergo self-organization. The possible self-organized states may or may not be spherically 
symmetric, although the velocity will be spherical symmetric. Self-organization implies that the plasma global 
conservation properties are sufficiently dominant to drive the system through dissipation, independently of the 
initial conditions. To find the possible self-organized states, we solve the ideal MHD equations, with 7 = 4/3, for 
self-similar solutions in time, where the temporal and spatial dependences of each physical variable are in sepa- 
rable form. We identify these self-similar solutions of ideal MHD to be the self-organized configurations reached 
by dissipation. We would like to emphasize the word ideal, which means that this process occurs without dissi- 
pation. Although self-organized states are physically generated by dissipation, they are part of the ideal MHD 
solutions. Dissipation provides the means by which an ideal state, of a particular energy density, is topologically 
transformed to another ideal state, of lower energy density. Using this time self-similar formulation, we can 
derive possible end states without time iteration, assuming a specified initial configuration. By circumventing 
intermediate stages, this self-similar approach, which is assumed to be independent of self-organization physical 
arguments, appears to be developing the required solutions for the problems waiting to be solved. 
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Under axisymmetry, we have established two solutions for equatorial ejection. The first solution has a finite 
azimuthal magnetic field, and represents a plasma torus with closed toroidal field lines. Beyond the plasma 
torus, the magnetic field lines are open and the plasma is not confined. This plasma torus provides the principle 
structure for the interacting-wind model. The second solution has a null azimuthal field, and represents a 
bipolar planetary nebula. Whether these self-organized structures continue to expand indefinitely depends on 
the integration constant H of the evolution function, which could have an asymptotically bounded track for a 
stable torus, and an unbounded track for an ever expanding bipolar nebulae. Additional nebula features, such 
as filaments and jets, could be accounted for with a fully three-dimensional self-similar MHD model. 

We in fact believe that this description of self-organized plasmas, using self-similar solutions, could represent 
a fundamental process in astrophysical ejection events. Apart from the equatorial ejection solutions of plasma 
torus and bipolar structure of planetary nebulae discussed here, the highly collimated polar ejection solutions 
could be relevant to extragalactic jets, should we consider these objects to be ejection events (Tsui and Serbeto 
2007). The collimated polar ejection mechanism could also be relevant to some asymmetric supernovae, with 
recoil on the neutron star, and to quasar ejection models of active galactic nucleus in cosmology. 

Acknowledgements. The author is deeply grateful to Prof. Akira Hasegawa for the very essential concept of self-organization in 
fluids and plasmas. 
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Fig. 1. The function Rq(z) = A zji(z) of the homogeneous solution is plotted as a function of z — arj to show 
the radial domains of the axisymmetric plasma structures. 




Fig. 2. The power series solution i?i(i 2 ) (z) of the particular solution, with A = 1 x 10 7 , is plotted as a function 
of z to show the plasma pressure effect on the radial domains of the axisymmetric plasma structures. 




Fig. 3. The power series solution i? 1 ( 9 )(z) of the particular solution, with A = 1 x 10 7 and n 3 — 16, is plotted as 
a function of z to show the plasma pressure effect on the radial domains of the axisymmetric plasma structures. 




Fig. 4. The power series solution i? 1 ( 6 )(z) of the particular solution, with A = 1 x 10 7 and k 3 = 16, is plotted as 
a function of z to show the plasma pressure effect on the radial domains of the axisymmetric plasma structures. 




Fig. 5. The power series solution i? 1 ( 3 )(z) of the particular solution, with A = 1 x 10 7 and k 3 = 16, is plotted as 
a function of z to show the plasma pressure effect on the radial domains of the axisymmetric plasma structures. 
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Fig. 6. The power series solution R^^z) of the particular solution, with A = 1 x 10 7 and k 3 = 16, is plotted as 
a function of z to show the plasma pressure effect on the radial domains of the axisymmetric plasma structures. 
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Fig. 7. The magnetic field lines using the homogeneous solution Rq(z) are shown with C = +0.05, C = +0.20, 
C = +0.40, C = +0.60, C = +0.80, and C = +1.00 from the outer to inner contours respectively. The axes are 
labelled in x = cos 8, and the contribution of the plasma pressure represented by R\{z) is neglected here. 




Fig. 8. The asymptotically bounded evolution function, the lower curve, is plotted as a function of the normalized 
time r with \H\/a = 0.5, whereas the unbounded evolution function, the upper curve, is plotted with H/a = 0.1. 
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Fig. 9. The explicit radial dependence of pressure profile pi(r]) is shown as a function of the normalized radial 
label 77/77* through the thick line. The explicit radial dependence of mass density profile pi(rj) is also shown 
through the thin line for comparisons. 
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Fig. 10. The full, explicit and implicit, radial mass density profile p(rj, P) using the homogeneous solution Rq(z) 
is shown as a function of z. 




Fig. 11. The mass density contours p{r],P) with C = are shown with C = +200, C = +300, C = +400, and 
C = +500 from the outer to inner contours respectively. 
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Fig. 13. The full, explicit and implicit, radial mass density profile p(rj, P) is shown as a function of z. 




Fig. 14. The magnetic field lines are shown with C = +0.0, C = +0.1, C = +0.3, C = +0.5, and C = +0.7 
from the outer to inner contours respectively. 



